#include <fstream>
#include <stdlib.h>
#include <math.h>

#include "Matrix.hpp"

using std::ifstream;
using std::ostream;

Matrix::Matrix( int n, int m ) : _n( n ), _m( m )
{
  createMatrix( &_valuesGPU, _m, _n );
  _values = new float[ _n*_m ];
}

void Matrix::print( ostream& out )
{
  
  out << _n << ' ' << _m << '\n';
  
  for( int i=0; i<_n; i++ )
  {
    for( int j=0; j<_m; j++ )
      out << _values[ i*_m+j ] << ' ';
    out << '\n';
  }
  
}

void Matrix::clear()
{
  
  delete _values;
  releaseMatrix( _valuesGPU );
  
}

Matrix Matrix::generate( int size )
{
  
  Matrix lMatrix = generateL( size );
  Matrix uMatrix = generateU( size );
  
  uploadMatrix( lMatrix._valuesGPU, lMatrix._values, lMatrix._m, lMatrix._n );
  uploadMatrix( uMatrix._valuesGPU, uMatrix._values, uMatrix._m, uMatrix._n );
  Matrix ret = lMatrix.mul( uMatrix );
  downloadMatrix( ret._valuesGPU, ret._values, ret._m, ret._n );
  
  lMatrix.clear();
  uMatrix.clear();
  
  return ret;
  
}

Matrix Matrix::generateL( int size )
{
  
  Matrix ret( size, size );
  
  for( int i=0; i<size; i++ )
    for( int j=0; j<size; j++ )
    {
      if( i<j )
	ret._values[ i*size+j ] = 0;
      else if( i==j )
	ret._values[ i*size+j ] = 1;
      else
        ret._values[ i*size+j ] = static_cast< float >( rand() )/RAND_MAX*20 - 10;
    }
    
 return ret;
  
}

Matrix Matrix::generateU( int size )
{
  
  Matrix ret( size, size );
  
  for( int i=0; i<size; i++ )
    for( int j=0; j<size; j++ )
    {
      
      if( i<j )
        ret._values[ i*size+j ] = static_cast< float >( rand() )/RAND_MAX*20-10;
      else if( i==j )
      {
	float value = static_cast< float >( rand() )/RAND_MAX*20-10;
	if( fabs( value ) < 1e-2 )
	  value += 1e-1 * ( value>=0 ? 1 : -1 );
	
	ret._values[ i*size+j ] = value;
      }
      else
	ret._values[ i*size+j ] = 0;
      
    }
  
  return ret;
  
}

Matrix Matrix::generateIdentity( int size )
{
  
  Matrix ret( size, size );
  generateIdentityMatrix( ret._valuesGPU, size );
  return ret;
  
}

Matrix Matrix::load( const char* filename )
{
  
  ifstream in( filename );
  
  int n, m;
  
  in >> n >> m;
  Matrix ret( n, m );
  
  for( int i=0; i<n; i++ )
    for( int j=0; j<m; j++ )
      in >> ret._values[ i*m+j ];
    
  return ret;
  
}

Matrix Matrix::mul( Matrix matrix )
{
  
  Matrix ret( _n, matrix._m );
  mul( matrix, &ret );
  
  return ret;
  
}

void Matrix::mul( Matrix matrix, Matrix *ret )
{
  
  mulMatrix( ret->_valuesGPU, _valuesGPU, matrix._valuesGPU, _n, _m, matrix._m );
  
}

Matrix Matrix::copy()
{
  Matrix ret( _n, _m );
  copy( &ret );
  
  return ret;
}

void Matrix::copy( Matrix *result, int yOffset, int xOffset )
{
  
  copyMatrix( _valuesGPU, result->_valuesGPU, _m, _n, result->_m, result->_n, xOffset, yOffset );

}

Matrix Matrix::submatrix( int xOffset, int yOffset, int xSize, int ySize )
{
  
  Matrix ret( ySize, xSize );
  
  getSubmatrix( ret._valuesGPU, _valuesGPU, xOffset, yOffset, xSize, ySize, _m, _n );
  
 return ret;
    
}

int Matrix::N()
{
  
  return _n;
  
}

void Matrix::doLUStep( Matrix lMatrix, int baseLine )
{
  
  makeLUStep( _valuesGPU, lMatrix._valuesGPU, baseLine, _m, _n );
  return;
  
}

void Matrix::doLU( Matrix lMatrix )
{
  
  for( int i=0; i<_n-1; i++ )
    doLUStep( lMatrix, i );
  
}

const int blockSize = 50;

void Matrix::doLUBlock( Matrix lMatrix )
{
  
  upload();
  lMatrix.upload();
  
  for( int i=0; i<_n; i+= blockSize )
   doLUBlockStep( lMatrix, i );
  
  download();
  lMatrix.download();
  
}

void Matrix::doLUBlockStep( Matrix lMatrix, int offset )
{
  
  Matrix subU(1,1), subL(1,1);
  
  if( offset+blockSize <= _n )
  {
    subU = submatrix( offset, offset, blockSize, blockSize );
    subL = generateIdentity( blockSize );
  }
  else
  {
    subU = submatrix( offset, offset, _n-offset, _n-offset );
    subL = generateIdentity( _n-offset );
  }
  subU.doLU( subL );
  
  subU.copy( this, offset, offset );
  subL.copy( &lMatrix, offset, offset );
  
  if( offset+blockSize >= _n )
    return;
  
  Matrix u2 = submatrix( offset+blockSize, offset, _m-offset-blockSize, blockSize );
  u2.solveSystemRight( subL );
  u2.copy( this, offset, offset+blockSize );
  
  Matrix l2 = submatrix( offset, offset+blockSize, blockSize, _n-offset-blockSize );
  l2.solveSystemLeft( subU );
  l2.copy( &lMatrix, offset+blockSize, offset );
  
  Matrix lu12 = l2.mul( u2 );
  subMatrix( lu12._valuesGPU, _valuesGPU, lu12._m, lu12._n, _m, _n, offset+blockSize, offset+blockSize );
  
  subL.clear();
  subU.clear();
  l2.clear();
  u2.clear();
  lu12.clear();
  
  clearUp( _valuesGPU, _n, offset+blockSize, offset, offset+blockSize ); 
  
}

void Matrix::solveSystemRight( Matrix lMatrix )
{
  for( int k=0; k<_n-1; k++ )
    solveSystemRightStep( _valuesGPU, lMatrix._valuesGPU, _m, _n, lMatrix._m, k );
  
  return;
    
}

void Matrix::solveSystemLeft( Matrix uMatrix )
{
  
  for( int i=0; i<_m; i++ )
  {
    
    scaleRow( _valuesGPU, uMatrix._valuesGPU, _m, _n, uMatrix._m, uMatrix._n, i );
    
    solveSystemLeftStep( uMatrix, i );
  }
}

void Matrix::download()
{
  
  downloadMatrix( _valuesGPU, _values, _m, _n );
  
}
void Matrix::upload()
{
  
  uploadMatrix( _valuesGPU, _values, _m, _n );
  
}

void Matrix::solveSystemLeftStep( Matrix uMatrix, int offset )
{
  
  systemRightStep( _valuesGPU, uMatrix._valuesGPU, _m, _n, uMatrix._m, offset );

}

bool Matrix::isEqual( Matrix matrix )
{
  
  int count = 0;
  for( int i=0; i<_n;i++ )
    for( int j=0; j<_m;j++ )
    {
      float delta = fabs( _values[ i*_m + j ] - matrix._values[ i*_m + j ] );
      if( delta > 0.1 )
	count++;
    }
  
  if( count > _n*_m*0.05 )
  {
    std::cout << count << '\n';
    return false;
  }
  
  return true;
  
}